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1.  Introduction 


Quantitative  understanding  of  laser-tissue  interaction  requires  a  thorough  knowledge  of 
the  distribution  and  propagation  of  light  intensity  in  biological  media.  This  is  particularly 
important  when  detennining  the  distribution  of  light  within  ocular  tissue  and  specifically 
in  determining  the  irradiance  at  the  retina  to  a  beam  of  light  incident  on  the  eye.  The  light 
distribution  intensity  can  be  obtained  from  the  solution  of  the  Helmholtz  equation  using 
the  slowly-varying  envelope  formalism.1  The  beam-propagation  method  is  at  present  the 
most  widely  used  tool  employed  in  this  study.  The  paraxial  (Fresnel)  approximation  is 
especially  popular  largely  owing  to  its  numerical  speed  and  simplicity. 

However,  the  paraxial  approximation  severely  limits  the  formalism  in  the 
following  two  respects.  First,  beams  containing  appreciable  Fourier  components  at  angles 
of  more  than  a  few  degrees  from  the  propagation  axis  will  experience  substantial  phase 
errors;  thus  the  method  cannot  treat  wide-angle  propagation.  Second,  beams  propagating 
through  regions  with  indices  of  refraction  that  differ  by  more  than  a  few  percent  from  the 
input  reference  index  will  also  suffer  serious  phase  distortion. 1 

The  Pade  approximant  scheme,  which  includes  the  paraxial  approximation, 
eliminates  some  of  the  error  in  the  paraxial  beam  propagation  model.  The  Pade 
approximant  scheme  is  also  computationally  efficient,  as  compared  to  an  approach  based 
on  the  method  of  lines.1  In  this  report,  we  develop  an  application  of  the  finite-difference 
method  to  the  Pade  approximant  formalism  for  optical  beam  propagation.  This  approach 
offers  substantial  improvements  in  both  wide-angle  propagation  and  index  variation 
tolerance  while  incurring  (in  the  two-dimensional  case)  only  a  modest  numerical 
penalty.1  Also,  employing  a  non-linear  map  from  a  infinite  physical  space  to  an  finite 
computational  space  avoids  spurious  reflections  from  the  computational  window  edge 
and  improves  computational  efficiency. 

In  addition,  a  model  sufficient  enough  to  accurately  describe  laser  propagation  in 
biological  media  must  also  address  time-dependent  optical-thermal  effects  including 
thennal  lensing  as  well  as  linear  and  non-linear  optical  effects.  Thermal  lensing  is  the 
gradient  of  the  refractive  index  caused  by  nonunifonn  heating.  Linear  optical  effects 
include  a  linear  correlation  between  the  absorption  and  scattering  coefficients  and  the 
temperature  of  the  medium.  Non-linear  optical  effects  include  non-linear  absorption, 
which  captures  the  dependence  of  the  absorption  coefficient  on  the  intensity  of  the 
incident  light.  Motamedi,  et  al.2  and  Lin,  et  al.3  have  shown  that  biological  materials  (e.g. 
tissue  and  blood)  possess  strong  thermal  lensing  effects.  Walsh  and  Cummings4  report 
that  moderate  transient  changes  of  the  absorption  coefficient  are  observed  in  IR-laser 
ablation  at  wavelengths  around  3  pm  (i.e.  close  to  the  absorption  peak  of  tissue  water). 
Lin  et  al.5  suggest  that  the  measured  transmittance  and  reflectance  decrease  and  increase, 
respectively  as  the  temperature-dependent  scattering  coefficient  increases.  Vodopyanov6 
shows  that  the  rise  of  the  water  temperature  during  the  ablation  process  can  induce  a 
decrease  of  the  absorption  coefficient  by  as  much  as  one  order  of  magnitude.  Stolarski, 
et  al.7  have  reported  non-linear  absorption  effects  of  melanin-doped  samples;  these 
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studies  are  of  particular  interest  to  laser  bioeffects  in  the  retina,  for  which  melanin  plays  a 
key  role  in  the  absorption  of  laser  energy. 

Moreover,  these  effects  are  interrelated.  A  sample  illuminated  with  light  will 
absorb  some  of  the  energy  as  quantified  by  the  local  intensity  of  light  and  local 
absorption  coefficient,  which  in  turn  depends  on  the  local  temperature  and  the  local 
recent  history  of  the  intensity.  The  absorbed  energy  causes  a  temperature  rise  that  induces 
a  thermal  lensing  effect.  This  thermal  lensing  effect  along  with  the  temperature- 
dependent  scattering  coefficient  changes  the  distribution  of  the  light.  Also,  Vodopyanov6 
shows  that  the  increase  in  temperature  of  water,  a  large  component  of  most  tissues,  is 
associated  with  a  shift  of  the  absorption  peak  towards  shorter  wavelengths.  These 
interdependencies  of  temperature  and  light  distribution,  though  seemingly  complex,  are 
capable  of  being  modeled  by  a  numerical  beam-propagation  model  coupled  to  a  time- 
dependent  optical-thermal  model.  Both  models  can  be  developed  independently  and 
coupled  in  a  time-slicing  scheme.  Moreover,  a  time-independent  beam-propagation 
model  depends  only  on  the  spatial  refractive  index  of  the  sample.  Therefore,  the  specific 
aim  of  this  report  was  to  develop  a  time-independent  finite-difference  model  of  light 
distribution  in  cylindrical  geometry  that  can  be  coupled  to  a  time-dependent  optical- 
thermal  model. 

As  an  application  of  the  above  mentioned  model,  the  irradiance  at  the  retina  from 
laser  light  entering  the  eye  can  be  predicted.  Upon  verification  of  the  model,  possibly 
with  a  z-scan,  the  predicted  incidence  light  at  the  retina  could  be  used  to  establish  new 
maximum  permissible  exposure  (MPE)  limits. 


2.  Non-Paraxial  Scalar  Wave  Equation 


Consider  the  scalar  Helmholtz  equation  in  cylindrical  geometry  obtained  by  using  the 
slowly-varying  envelope  formalism: 


a 

dH(r,z) 

dr 

-  dr 

dz 


+ 


d2H(r2Z)  +  k2Q  [ i2( r)  -  n2  }/(r,zJ  =  0 


(1) 


where  k  =  k0n0  with  k0  the  free-space  propagation  constant,  n0  the  (input)  reference 
refractive  index.  We  map  the  radial  coordinate  r  of  infinite  support  [0,oo]  to  a  variable  of 
finite  support  [0,zr/  2]  through  the  mapping  function 

r  =  ij  tan  p  (2) 

Where  77  is  the  coordinate  scaling  parameter  that  allows  us  to  control  the  sampling.  Such 
a  representation  maps  onto  the  entire  physical  space,  eliminating  boundary  conditions 
that  could  lead  to  spurious  reflections,  as  well  creating  a  unifonn  grid  from  a  non- 
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uniform  grid  such  that  a  finer  grid  spacing  is  placed  around  the  origin.8,9  The  latter  is 
quite  useful  for  modeling  beam  propagation  with  cylindrical  symmetry  about  the  origin. 


Substitution  of  Eq.  (2)  into  Eq.  (1)  yields 

dH i  82H  _  iP 

dz  2k  dz 2  2k 

where  the  operator  P  is  defined  by: 

P  =  k2  [n 2  (r  )  -  n]  ]  +  C,°S  P  (2  cos 2  p  - 1)  ^  . 

77“  sin/}  op 


Formally  rewriting  Eq.  (3)  in  the  form: 


a H 

dz 


iP 


2k 


2k  dz 


H 


suggests  the  recursion: 


d_ 

dz 


H  = 


iP 

2k 


1  c 
2k  dz 


-H 


1- 


(3) 


(4) 


(5) 


(6) 


With 


=  0  ,  we  arrive  at  the  expression: 


dH 

Ik 


iP  |  iP2 
2k  +  4 kz  H 

,  3  P  P2 

1  H - v  ^ - Y 

4  k2  16  k4 


(V) 
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Equation  (7)  is  the  (2,2)  Pade  approximant  for  the  Helmholtz  operator.  Hadley1  shows 
the  relative  phase  error  incurred  by  the  use  of  the  (2,2)  Pade  approximate  operator  is  less 

than  0.2  for  up  to  65°  for  an  incident  plane  wave. 


3.  Finite-Difference  Model 

Discretization 

Finite-difference  equations  may  be  derived  from  Eq.  (7)  yielding: 


iAz 


D(Hj  -  H™  )  =  —  N(HJ 


+  ET+1) 


(8) 


where  N  and  D  are  the  numerator  and  denominator,  respectively,  of  Eq.  (7).  The 
superscript  in  Eq.  (8)  designates  the  z  position  while  the  subscript  designates  the  radial 
position.  Substituting  for  N  and  D  for  the  numerator  and  denominator  in  Eq.  (8), 
respectively,  from  Eq.  (7)  yields: 


H"/+l  +  aPH'"+l  +  PP2H ;+1  =  H™  +  yPHJ  +  SP2H"‘ 

where 


3  iAz  1  iAz 

~4k2~^k  6k*  8 ¥ 


(9) 


(10) 


The  use  of  centered  spatial  differencing  results  in  the  following  form  for  the  operator  P 
for  the  axially  symmetric  case: 


PH™ 


A>H'U 


+  BJHJ  +  CH 


m 

j+ 1 


and 


P2HJ 


Aj  +  B™_XH™_X  +  Cj_xHJ  ) 

+  BJ(AjH'^+BJHJ  +CjHJ+l) 

+  Cj(Aj+lHj  +  BJ+1H™+1  +  Cj+1HJ+2) 


(ID 


(12) 


where 


cos4  pj  l 

r  (A pf 


cos3  pj  1 
77 2  sin  Pj  2(A p)2 


(2  cos2  Pj  -1) 
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(13) 


B"j  =  tf\n2(pj,zm)-nl\-2 


cos  p.  l 


r 


M2 


cos  p  1  cos  p  1  /  2  \ 

c,  = - —7 — 77+  ,  .  • - 7 — 7712COS  p.  -l)  . 

7  (A/?)  7‘sinp  ■  2(A/?) 


Substituting  Eqs.  (11)  and  (12)  into  Eq.  (9)  yields 

PAjA^H”*  +{aAj  +  +  pAj  B"(+')h  n(+( 

+  (l  +  a5;+1  +  +  j3AJ+lCj  )h’J+' 

HaCj  +  fiBJ+'C  j  +/3B$CJ)H$  +flCJCJ,lH% 

=  rAjA^H^  +  +  &i(s;  Iff-, 

+  (l  +  yBn;  +  <54  .Cy_!  +  SB1"  BJ  +  8Aj+xC}  )hJ 
+  {>C,  +SB-C ,  +«-,C;)ff-|  +XJC/tlH-’,1 

It  is  clear  from  Eq.  (14)  that  the  system  is  pentadiagonal. 


(14) 


4.  Boundary  Conditions 

Assuming  axial  symmetry  in  the  radial  direction  (i.e.  p  f  =  -pj  and  H”,  =  HJ)  we 
arrive  at  the  following  boundary  conditions  for  j  =  0  and  j  =  1: 


1  +  aB(,+1  +pB"'+1B(,+l  + 


+ 


+ 


2p  cos4  px 
7  4  (A p)4  r/4(Apf 

^  2a  2p  x  2(3  ,i+1A 

H 77 - 77“  On  - 77 - 77  B 


(3  cos3  /?j(2cos2  px -1) 


sin  px 


H, 


m+ 1 


7 


!(A pf  72(A/i) 


2  ^0 


7 


!(Ap) 


2  ~  1 


//, 


*2+1 


J 


2 [3  [3  cos3  A  (2 cos2  A -1) 

—  COS  A  H 


74(Ap) 

i +rBo  +sb"b;  + 


sin  px 


7 4  (A py 
28  cos4  a 
7 4  (A/i  )4  74(Aa)3 


H: 


t+1 


/ 


A  cos3  a  (2  cos 1  a  ~  1) 


sin  a 


///; 


2y 


-  +  - 


2A 


-s;+- 


2A 


A 


+ 


7 2  (A/?)2  72(A/?)2  °  72(Ap)  y 

28  4  A  cos3  a  (2  cos2  a  - 1) 

?M^C0S  Pl+VM 


-B( 


H 


sin  px 


j 


(15) 
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and 


(oA,  +  /3AXB;+X  + 1 3AxB[n+x)H , 


m+ 1 
0 


I  +  aB™+l  +  J3A 


A  4 

2  cos  px 


+  pB™+lB™+l  +  /U2C, 


r(/±p) 

+  (aci  +  +  fjB'"+ict  )//2ml  +  pcfji 

=  (/Al  +SA1B’n  +SA1B”)hi 

9  COS4  /9  ^ 

+  1  +  yB!n  +  M  — — ^  +  SB!nB!n  +  &LC, 

+  (yC,  +  (ffif'C,  +  6B'X\  )//2'  +  SCxC2H'" 


H 


m+l 


m+l 

3 


ii : 


(16) 


5.  Computer  Program 

The  finite-difference  model  described  above  was  implemented  in  C++  and  is 
included  in  Appendix  A.  The  pentadiagonal  solver  is  based  on  LU  decomposition  as 
described  in  Ref.  10.  At  the  time  this  report  was  written,  the  code  would  compile  but  not 
yet  run  as  expected.  This  is  most  likely  due  to  a  logical  error  in  the  code  which  should  be 
a  minor  fix. 


6.  Conclusions 

In  the  report,  a  finite-difference  model  for  beam  propagation  within  cylindrical 
geometry  is  developed  using  a  slowly  varying  envelope  fonnalism  for  the  scalar 
Helmholtz  equation.  A  finite-difference  computer  model  for  the  (2,2)  Pade  approximant 
to  the  Helmholtz  operator  was  used  as  it  is  capable  of  handling  wide-angle  propagation 
and  refractive  index  variation,  while  still  maintaining  numerical  speed  and  simplicity  A 
transfonnation  from  the  infinite  physical  space  to  a  finite  computational  space  was  also 
employed,  eliminating  the  boundary  conditions  at  the  edge  of  the  computational  window 
and  possible  spurious  reflections. 

This  computer  model,  developed  in  C++,  is  able  to  couple  to  a  thennal  model  in  a 
time-slicing  scheme  and  could  be  used  to  model  the  time-dependent  light  distribution, 
including  thermal  lensing  effects,  within  the  eye  as  well  as  the  irradiance  incident  at  the 
retina.  Upon  verification  of  the  model  with  experiment,  for  example  by  a  z-scan, 
predicted  light  distribution  values  within  the  eye  could  be  used  for  practical  applications. 
In  particular,  the  predicted  irradiance  at  the  retina  of  laser  light  incident  on  the  human  eye 
could  be  used  to  establish  new  MPE  limits. 
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Appendix  A 


#inc  lude<stdio .  h> 

#inc  lude<stdlib .  h> 

#inc  lude<iostream.  h> 

#include<complex.h> 

#inc  lude<math .  h> 

#include<fstream.h> 

double  get_n(complex<double>  r,  complex<double>  z); 

void  get_EO(complex<double>  nnax,  complex<double>  delr,  int  rn,  complex<double>  soln); 

void  pentadiagonalsolver(int  n  ,  complex<double>  g,  complex<double> 
h,complex<double>  d,  complex<double>  e,  complex<double>  ,  complex<double> 
b,  complex<double>  soln); 


double  get_n(complex<double>  r,  complex<double>  z) 

{ 


complex<double>  deln; 
double  n; 
deln  =  0.00005; 
double  lens  =  0.0; 

n  =  ((1.5  +  1.0*(-lens))  +  (lens  /(1+  real((r)/(del_n))*real((r)/(del_n))))); 
return(n); 


} 
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void  get_EO(complex<double>  nnax,  complex<double>  delr,  int  rn,  complex<double>  soln) 

{ 

int  j; 

for(j=0;j<  m;j++) 

{ 

soln[j]  =  1  *  (std::exp(-(pow(((  (j)*  real(delr))  /  real((rmax)*0.1)  ),2)))  ); 

} 


} 

void  pentadiagonalsolver(int  n,  complex<double>  g,  complex<double>  h, 
complex<double>  d,  complex<double>  e,  complex<double>  f,  complex<double>  b, 
complex<double>  soln) 

{ 


intk; 

complex<double>  alpha[n]; 
complex<double>  gam[n-l]; 
complex<double>  delta[n-2]; 
complex<double>  bet[n]; 
complex<double>  c[n]; 
alpha[0]  =  d[0]; 
gam[0]  =  e[0]/alpha[0]; 
delta [0]  =  f[0]  /  alpha[0]; 
bet[l]  =  h[l]; 

alpha[l]  =  d[l]  -  (bet[  1  ]  *  gam[0]); 
gam[l]  =  (e[l]  -  (bet[l]  *  delta  [  0] )  )/alpha  [  1  ] ; 
delta[l]  =  f[  1  ]/  alpha[  1  ] ; 
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for(k=2  ;k<=n-3  ;k++) 


bet[k]  =  h[k]  -  (  g[k]*gam[k-2]); 

alpha[k]  =  d[k]  -  (g[k]  *  delta[k-2])  -  (bet[k]  *  gam[k-l]); 
gam[k]  =  (  e[k]  -  (bet[k]*  delta [k-l]))/alpha[k]; 
delta[k]  =  f[k]/alpha[k]; 


} 

bet[n-2]  =  h[n-2]  -  (g[n-2]  *  gam[n-4]); 

alpha[n-2]  =  d[n-2]  -  (g[n-2]*delta[n-4])  -  (bet[n-2]  *  gam[n-3]); 
gam[n-2]  =  (  e[n-2]  -  (bet[n-2]  *  delta[n-3])  )  /  alpha[n-2]; 
bet[n-l]  =  h[n-l]  -  (g[n-l]  *  gam[n-3]); 

alpha[n-l]  =  d[n-l]  -  (g[n-l]  *  delta[n-3])  -  (bet[n-l]  *  gam[n-2]); 

c[0]  =  b[0]  /  alpha[0]; 

c[l]  =  (b[  1  ]  -  (bet[  1  ] * c [0] ))/ alpha [  1  ] ; 

for(k=2  ;k<n;k++) 

{ 

c[k]  =  (b[k]  -  (g[k]*c[k-2])  -  (bet[k]*c[k-l]))/alpha[k]; 

} 

soln[n-l]  =  c[n-l]; 

soln[n-2]  =  c[n-2]  -  gam[n-2]  *soln[n-l]; 


for(k=n-3;k>=0;k— ) 
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{ 

soln[k]  =  c[k]  -  (gam[k]  *  soln[k+l])  -  (delta[k]  *  soln[k+2]); 

} 


int  main() 

{ 

complex<double>  lambda  =  0.000001; 
complex<double>  kO  =  ((  2  *  3.14159  )/  lambda); 
complex<double>  nbar  =1.5; 
complex<double>  k  =  kO  *  nbar; 
complex<double>  rhomax  =  0.00001; 
complex<double>  zmax  =  0.00004; 
complex<double>  delz  =  0.0000001; 
complex<double>  delrho  =  0.0000001; 
complex<double>  rhonl; 

rhonl  =  rhomax  /  delrho; 
int  rhon  =  real(rhonl)  +  1 ; 
complex<double>  i; 
imag(i)  =  1; 
real(i)  =  0; 

complex<double>  znl; 

znl  =  zmax  /  delz; 

I  nt  zn  =  real(znl)  +  1; 
complex<double>  eta  =  1.0; 
complex<double>  DLa[rhon]; 
complex<double>  DLb[rhon]; 
complex<double>  DLc[rhon]; 
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complex<double>  DLd[rhon]; 
complex<double>  DLe[rhon]; 
complex<double>  DRa[rhon]; 
complex<double>  DRb[rhon]; 
complex<double>  DRc[rhon]; 

complex<double>  DRd[rhon]; 
complex<double>  DRe[rhon]; 
complex<double>  soln[rhon]; 

intp; 
intj; 
int  m; 

complex<double>  A[3]; 
complex<double>  B[3]; 
complex<double>  C[3]; 
complex<double>  rho[rhon+l]; 
complex<double>  z[rhon]; 

rho[0]  =  0; 
z[0]  =  0; 

for(p=l  ;p<rhon;p++) 

{ 

rho[p]  =  rho[p-l]  +  delrho; 

} 

for(p=  1  ;p<zn;p++) 
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{ 


z[p]  =  z[P‘l]  +  delz; 


} 

rho[rhon]  =  rhomax  +  delrho;  //  Extra  Point  Needed  in  Loop  when  j=rhon-l 


complex<double>  alpha; 
complex<double>  beta; 
complex<double>  gama; 
complex<double>  delta; 


alpha  =  ( (3. 0/(4. 0  *  k  *  k))  -  ( (i  *  delz)/(4.0  *  k)  )  ); 

beta  =  ((1.0  /  (  16.0  *  k  *  k  *  k  *  k))  -  ((i  *  delz  )/  (8.0  *  k  *  k  *  k))  ); 

gama  =  ( (3. 0/(4. 0  *  k  *  k))  +  ( (i  *  delz)/(4.0  *  k)  )  ); 

delta  =  ((1.0  /  (  16.0  *  k  *  k  *  k  *  k))  +  ((i  *  delz  )/  (8.0  *  k  *  k  *  k))  ); 

double  n; 


get_E0(rhomax,  delrho,  rhon,soln); 

complex<double>  *01dE; 

OldE  =  soln; 

double  output  1,  output2; 
ofstream  outputfile; 
outputfile .  op  en( "  output .  txt" ) ; 
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for(j =0  ;j  <rhon;j ++) 

{ 

outputl  =  real(01dE[j]); 
output2  =  imag(01dE[j]); 
outputfde«output  1  <<"\t\t"«output2«"\n" ; 


} 

for(m=  1  ;m<zn;m++) 

{ 

//  Boundary  Conditions 

DLa[0]  =  0; 

DLb[0]  =  0; 
n  =  get_n(rho[0],  z[m]); 

B[0]  =  ( ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[0]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[l],  z[m]); 

B[l]  =  (  ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

DLc[0]  =  (1.0  +  (alpha  *  B[0])+  (beta  *  B[0]  *  B[0])+  ( (2.0  *  beta  *  (pow  ( (  cos  (real(rho[l]))),  4)))/((  pow  (eta,  4))*(  pow  ( 
delrho  ,  4  ))  )  )  -  ( (beta  *  (pow(  (cos  (real(rho[l]))),  3)  )  *  (2.0  *  (  pow  ((cos(real(rho[l])))  ,3  ))  )  -  1.0  )/  ( (pow(  eta,  4))  *  (  pow  ( 
delrho  ,  3))  *  sin(real(rho[l]))  )  )  ); 

DLd[0]  =  ( ((2.0  *  alpha)/(eta  *  eta  *  delrho  *  delrho))+  ((2.0  *  beta  *  B[0])/(eta  *  eta  *  delrho  *  delrho))+  ((2.0  *  beta  * 
B[l])/(eta  *  eta  *  delrho  *  delrho)) ); 
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DLe[0]  =  ( ( (2.0  *  beta  *  (  pow  ((cos  (real(rho[l]))),4  )))/(  pow(eta,  4)*  pow(delrho,4)  ))  +  ( (beta  *  (pow(  (cos(real(rho[l]))), 
3)  )  *  (2.0  *  (  pow  ((cos(real(rho[l])))  ,3  ))  )  -  1.0  )/  ( (pow(  eta,  4))  *  (  pow  (  delrho  ,  3))  *  sin(real(rho[l]))  )  )  ); 

DRa[0]  =  0; 

DRb[0]  =  0; 

n  =  get_n(rho[0],  z[m-l]); 

B[0]  =  ( ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[0]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

n  =  get_n(rho[l],  z[m-l]); 

B[l]  =  (  ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos(real(rho[l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

DLc[0]  =  (1.0  +  (gama  *  B[0])+  (delta  *  B[0]  *  B[0])+  ( (2.0  *  beta  *  (pow  ( (  cos  (real(rho[l]))),  4)))/((  pow  (eta,  4))*(  pow  ( 
delrho  ,  4  ))  )  )  -  ( (delta  *  (pow(  (cos  (real(rho[l]))),  3)  )  *  (2.0  *  (  pow  ((cos(real(rho[l])))  ,3  ))  )  -  1.0  )/  ( (pow(  eta,  4))  *  (  pow  ( 
delrho  ,  3))  *  sin(real(rho[l]))  )  )  ); 

DLd[0]  =  ( ((2.0  *  gama)/(eta  *  eta  *  delrho  *  delrho))+  ((2.0  *  delta  *  B[0])/(eta  *  eta  *  delrho  *  delrho))+  ((2.0  *  delta  * 
B[l])/(eta  *  eta  *  delrho  *  delrho))); 

DLe[0]  =  ( ( (2.0  *  delta  *  (  pow  ((cos  (real(rho[l]))),4  )))/(  pow(eta,  4)*  pow(delrho,4)  ))  +  ( (delta  *  (pow(  (cos 
(real(rho[l]))),  3)  )  *  (2.0  *  (  pow  ((cos(real(rho[l])))  ,3  ))  )  -  1.0  )/  ( (pow(  eta,  4))  *  (  pow  (  delrho  ,  3))  *  sin(real(rho[l]))  )  )  ); 

DLa[l]  =  0; 

A[0]  =  ( ( (pow(  (cos  (real(rho[0]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[0]))),  3  )  )  /  (eta  * 
eta  *  sin(  real(rho[0])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[0])))  ,2)  )  )  -  1.0  )  )  ; 
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A[l]  =  (  ( (pow(  (cos  (real(rho[l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[l]))),  3  )  )  /  (eta  * 

eta  *  sin(  real(rho[l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[l])))  ,2)  )  )  -  1.0  ))  ; 

A[2]  =  ( ( (pow(  (cos  (real(rho[2]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[2]))),  3  )  )  /  (eta  * 
eta  *  sin(  real(rho[2])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[2])))  ,2)  )  )  -  1.0  ))  ; 

n  =  get_n(rho[0],  z[m]); 

B[0]  =  ( ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[0]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

n  =  get_n(rho[l],  z[m]); 

B[l]  =  (  ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

n  =  get_n(rho[2],  z[m]); 

B[2]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[2]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

C[0]  =  ( ( (pow(  (cos  (real(rho[0]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((eos(real(rho[0]))),  3  )  )  /  (eta  * 

eta  *  sin(real(  rho[0])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos(real  (rho[0])))  ,2)  )  )  -  1.0  )  )  ; 

C[l]  =  (  ( (pow(  (cos  (real(rho[l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((cos(real(rho[l]))),  3  )  )  /  (eta  * 

eta  *  sin(real(  rho[l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[l])))  ,2)  )  )  -  1.0  )  )  ; 

C[2]  =  ( ( (pow(  (cos  (real(rho[2]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((cos(real(rho[2]))),  3  )  )  /  (eta  * 

eta  *  sin(real(  rho[2])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[2])))  ,2)  )  )  -  1.0  )  )  ; 

DLb[l]  =  ((A[l]  *  alpha)  +  (beta  *  A[l]  *  B[0])  +  (beta  *  B[l]  *  A[l])  ); 

DLc[l]  =  (  1.0  +  (B[  1]  *  alpha)  +  ( (  beta  *  A[l]  *  2.0  *  pow(  (cos(real(rho[l]))),4))/(  eta  *  eta  *  delrho  *  delrho)  )  +  (beta  * 

B[l]  *  B[l])  +  (  beta  *  C[l]  *  A[2])  ); 
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DLd[l]  =  ((C[l]  *  alpha)+  (beta  *  B[l]  *  C[l])+  (beta  *  C[l]  *  B[2])  ); 

DLe[l]  =  (beta  *  C[l]  *  C[2]); 
n  =  get_n(rho[0],  z[m-l]); 

B[0]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[0]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[l],  z[m-l]); 

B[l]  =  (  ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[2],  z[m-l]); 

B[2]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[2]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

DRa[l]  =  0; 

DRb[l]  =  ((A[l]  *  gama)  +  (delta  *  A[l]  *  B[0])  +  (delta  *  B[l]  *  A[l])  ); 

DLc[l]  =  (  1.0  +  (B[  1]  *  gama)  +  ( (  delta  *  A[l]  *  2.0  *  pow(  (cos(real(rho[l]))),4))/(  eta  *  eta  *  delrho  *  delrho)  )  +  (delta  * 
B[l]  *  B[l])  +  (  delta  *  C[l]  *  A[2]) ); 

DRd[l]  =  ((C[l]  *  gama)  +  (delta  *  B[l]  *  C[l])  +  (delta  *  C[l]  *  B[2])  ); 

DRe[l]  =  (delta  *  C[l]  *  C[2]); 

for(j  =2  ;j  <rhon;j  ++) 

{ 
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A[0]  =  (  (  (pow(  (cos  (real(rho[j-l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[j-l]))),  3  ) 
)  /  (eta  *  eta  *  sin(real(  rho[j-l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j-l])))  ,2)  )  )  -  1.0  )  )  ; 

A[l]  =  (  (  (pow(  (cos  (real(rho[j]  ))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[j]  ))),  3  )  )  / 

(eta  *  eta  *  sin(real(  rho[j]  )))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j]  )))  ,2)  )  )  -  1.0  )  ); 

A[2]  =  ( ( (pow(  (cos  (real(rho[j+l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))-(((  pow((cos(real(rho[j+l]))),  3 
)  )  /  (eta  *  eta  *  sin(real(  rho[j+l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j+l])))  ,2)  )  )  -  1.0  )  )  ; 

n  =  get_n(rho[j-l],  z[m]); 

B[0]  =  ( ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[j-l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[j],  z[m]); 

B[  1]  =  ( ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[j]  ))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[j+l],  z[m]); 

B[2]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[j+l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 

C[0]  =  ( ( (pow(  (cos  (real(rho[j-l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((cos(real(rho[j-l]))),  3  )  )  /  (eta 

*  eta  *  sin(real(  rho[j-l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j-l])))  ,2)  )  )  -  1.0  )  )  ; 

C[l]  =  ( ( (pow(  (cos  (real(rho[j]  ))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((cos(real(rho[j]  ))),  3  )  )  /  (eta  * 
eta  *  sin(real(  rho [j ]  )))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j]  )))  ,2)  )  )  -  1.0  )  )  ; 

C[2]  =  ( ( (pow(  (cos  (real(rho[j+l]))  ),  4))/  (eta  *  eta)  )  *  (  1.0/  (delrho  *  delrho)  ))  +  (((  pow((cos(real(rho[j+l]))),  3  )  )  / 
(eta  *  eta  *  sin(real(  rho[j+l])))  )  *  (1.0/  (2.0  *  delrho)  )  *  ( (2.0  *  (pow(  (cos  (real(rho[j+l])))  ,2)  )  )  -  1.0  )  )  ; 

DLa[j]  =  (beta  *  A[l]  *  A[0]); 
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DLb[j]  =  ((A[  1  ]  *  alpha)  +  (beta  *  A[l]  *  B[0])  +  (beta  *  B[l]  *  A[l]) ); 

DLc[j]  =  (  1.0  +  (B[  1  ]  *  alpha)  +  (beta  *  A[l]  *  C[0])  +  (beta  *  B[l]  *  B[l])  +  (beta  *  C[l]  *  A[2])); 

DLd[j]  =  ((C[l]  *  alpha)+  (beta  *  B[l]  *  C[l])+  (beta  *  C[l]  *  B[2]) ); 

DLe[j]  =  (beta  *  C[l]  *  C[2]); 
n  =  get_n(rho[j-l],  z[m-l]); 

B[0]  =(  ((k0  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos  (real(rho[j-l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[j],z[ni-l]); 

B[  1]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos(real(rho[j]  ))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
n  =  get_n(rho[j+l],  z[m-l]); 

B[2]  =  ( ((kO  *  kO)  *  ((n  *  n  )-(nbar  *  nbar)  ))  -  ( (  2.0  *  pow(  (cos(real(rho[j+l]))),  4))/(eta  *  eta  *  delrho  *  delrho  )  )); 
DRa[j]  =  0; 

DRb[j]  =  ((A[l]  *  gama)  +  (delta  *  A[l]  *  B[0])  +  (delta  *  B[l]  *  A[l])  ); 

DRc[j]  =  ( 1.0  +  (B[  1  ]  *  gama)  +  (delta  *  A[l]  *  C[0])  +  (delta  *  B[l]  *  B[l])  +  (delta  *  C[l]  *  A[2])); 

DRd[j]  =  ((C[l]  *  gama)  +  (delta  *  B[l]  *  C[l])  +  (delta  *  C[l]  *  B[2])  ); 

DRe[j]  =  (delta  *  C[l]  *  C[2]); 
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} 


DLe[rhon-2]  =  0; 

DLd[rhon-l]  =  0; 

DLe[rhon-l]  =  0; 

DRd[rhon-l]  =  0; 

DRe[rhon-l]  =  0; 

DRe[rhon-2]  =  0; 

//  Compute  RHS  Column  Vector 

complex<double>  D[rhon]; 

D[0]  =  (  DRc[0]  *  OldE[0] )  +  (  DRd[0]  *  01dE[l] )  +  (  DRe[0]  *  01dE[2]); 

D[l]  =  (  DRb[0]  *  OldE[0]  )  +  (  DRc[0]  *  01dE[l]  )  +  (  DRd[0]  *  01dE[2])  +  (  DRe[0]  *  01dE[3]  ); 


for(j =2  ;j  <rhon-2  ;j ++) 

{ 


DU]  =  (  DRaU]  *  OldEU-2]  )  +  (  DRbU]  *  01dE[j-l]  )  +  (  DRcU]  *  OldEU]  )+  (  DRdU]  *  OldEU+1]  )  +  (  DReU]  *  OldEU+2] 

); 


} 

D[rhon-l]  =  (  DRa[rhon-l]  *  01dE[rhon-3]  )  +  (  DRb[rhon-l]  *  01dE[rhon-2]  )  +  ( 

DRc[rhon-l]  *  01dE[rhon-l]  ); 

D[rhon-2]  =  (  DRa[rhon-2]  *  01dE[rhon-4]  )  +  (  DRb[rhon-2]  *  01dE[rhon-3]  )  +  (  DRc[rhon-2]  *  01dE[rhon-2]  )  +  (  DRd[rhon-2]  * 
01dE[rhon-l]  ); 
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//  Solve  System 

pentadiagonalsolver(rhon,  DLa,  DLb,  DLc,  DLd,  DLe,  D,  soln); 


OldE  =  soln; 

//  Write  to  Output  File 

for(j=0;j<rhon;j++) 

{ 

outputl  =  real(01dE[j]); 
output2  =  imag(01dE[j]); 
outputfde«output  1  <<"\t\t"«output2«"\n" ; 

} 

} 

outputfile.close(); 

return(O); 

} 
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